Spatio-temporal wave propagation in photonic crystals: a Wannier-function analysis 
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A general analysis of undistorted propagation of localized wavepackets in photonic crystals based on 
a Wannier-function expansion technique is presented. Different kinds of propagating and stationary 
spatio-temporal localized waves are found from an asymptotic analysis of the Wannier function 
envelope equation. 

PACS numbers: 41.20.Jb, 42.70.Qs, 42.25.Bs 



I. INTRODUCTION 



Spatio-temporal broadening of localized wavepackets 
with finite energy due to the effects of diffraction and 
dispersion is a universal and challenging phenomenon in 
any physical context involving wave propagation. If the 
finite energy constraint is left, special spatio-temporal 
waves with a certain degree of localization in space 
and/or in time, capable of propagating free of diffrac- 
tion and/or temporal dispersion, can be constructed. Lo- 
calized waves of this type include, among others, Bessel 
beams, focus-wave modes, X-type waves, pulsed Bessel 
beams 0, 0, S 0] ■ Though these waves can be only 
approximately realized in practice, several experiments in 
acoustic and optical fields have been reported so far show- 
ing nearly-undistorted localized wave propagation. As 
the existence of undistorted progressive localized waves 
in vacuum has been known since many years and lead 
to long-standing studies P, Q , with special attention de- 
voted toward their superluminal or subluminal character 
and to their finite-energy realizations, in the past few 
years these studies have been extended to dispersive op- 
tical media [1,^1^ , and remarkably the spontaneous gen- 
eration of localized and nonspreading wavepackets medi- 
ated by optical nonlinearities has been predicted 0] and 
experimentally observed iZ] using standard femtosecond 
pulsed lasers. Very recently, in a few works H, 0, 
the issue of spatial or spatio-temporal wave localization 
in periodic media has been addressed, and the possibility 
of exploiting well-established anomalous diffractive and 
dispersive properties of photonic crystals (PCs) 0, 0| 
to induce novel spatio-temporal wave localization mecha- 
nisms has been proposed. Specifically, these studies have 
been concerned with localization of Bose-Einstein con- 
densates in a one-dimensional optical lattice without any 
trapping potential 8], with two-dimensional (2D) spa- 
tial Bessel X waves in weakly-coupled 2D waveguide ar- 
rays showing bi-dispersive properties and with three- 
dimensional (3D) out-of-plane X-wave localization in 2D 
PCs 10]. Spatio-temporal waves considered in these 
works rely on some specific models and often use ad-hoc 
approximations, e.g. reduced coupled-mode equations, 
paraxiality, weak-coupling limit, continuum approxima- 
tions. So far, a general framework to capture spatio- 
temporal wave localization and propagation in PCs and 



the derivation of a general wave equation, valid regard- 
less of the specific system under investigation and with a 
wide range of applicability, is still lacking. 
The aim of this work is to provide a general analytical 
framework to study spatio-temporal wave propagation in 
2D and 3D PCs based on the use of Wannier-functions, 
which have been introduced in the context of PCs to 
treat localized modes, such as the bound states of im- 
purities or lattice defects 0, 0| . A general asymptotic 
analysis of the envelope equation for the Wannier func- 
tions allows one to capture the existence and properties 
of localized nonspreading wavepackets in PCs in terms 
of localized solutions of canonical wave equations, such 
as the Schrodinger equation, the Helmholtz equation and 
the Klein-Gordon equation. 



II. WANNIER FUNCTION ENVELOPE 
EQUATION 

The starting point of the analysis is provided by the 
vectorial wave equation for the magnetic field H = 
H(r, t) in a PC with a periodic relative dielectric con- 
stant e(r). 
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where c is the speed of light in vacuum. To study the 
propagation of a spatio-temporal wavepacket, we can 
adopt the method of the Wannier function expansion, 
which is commonplace in the study of the quasi-classical 
electron dynamics in solids 15. 16j and recently ap- 
plied to study localized modes and defect structures in 
PCs with defects [Hill. We refer exphcitly to a 3D 
PC structure, however a similar analysis can be devel- 
oped for a 2D PC. Let us first consider the monochro- 
matic Bloch-type solutions to Eq.(l) at frequency w, 
H(r, t) = Hk.n(r) exp(— iwt), where k lies in the first 
Brillouin zone of the reciprocal k space, w = u;„(k) is 
the dispersion curve for the ri-th band, and Hk,ri(r) are 
the band modes, satisfying the condition Hk,,i(r + R) = 
Hk,n(r) exp(jk • R) for any lattice vector R of the peri- 
odic dielectric function. The Bloch functions Hk.n (r) are 
normalized such that (Hk'.n' |Hk.n) — Vbz i5n,n"5(k' — k), 
where Vbz — (27r)^/F is the volume of the first Brillouin 
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zone in the reciprocal space and V is the volume of the 
real-space unit cell. For each band of the PC, one can 
construct a Wannier function W„(r) as a localized su- 
perposition of Bloch functions of the band according to: 



W„(r) = [ dk Hk,„(r). 

Vbz Jbz 



(2) 



In the superposition, the phase of Bloch functions Hk,n 

can be chosen such that the Wannier function W„(r) 
is strongly localized around r = with an exponential 
decay far from r = 0. The Wannier functions satisfy the 
orthogonality conditions (W„'(r — R')|W„(r — R)) = 
Sn,n' Sr,r', and the following relationship can be easily 
proven: 



(W„.(r-R') 
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X^R,^n,Rexp(ik-R). We then look for a spatio-temporal 
wavepacket, which is a solution to Eq.(l), as a super- 
position of translated Wannier functions localized at the 
different lattic;e points R of the periodic structure, with 
amplitudes f(R,t) that depend on the lattice point R 
and can vary in time, i.e. we set: 



H(r,i)=^/(R,<)W4r-R). 



(4) 
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Note that, as we consider a pure periodic structure with- 
out defects and neglect perturbation terms in Eq.(l) (e.g. 
nonlincaritics), coupling among different bands does not 
occur and in Eq.(4) the sum can be taken over a sin- 
gle band, of index n. Coupled-mode equations for the 
temporal evolution of the amplitudes /(R, t) of Wannier 
functions at different lattice points can be obtained af- 
ter substitution of Eq.(4) into Eq.(l), taking the scalar 
product with W„(r — R) and using the orthogonality 
conditions of Wannier functions, together with Eq.(3). 
One obtains: 



(5) 
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The solution to the couplcd-modc equations (5) can be 
expressed as /(R, t) = /(r = R, t), where the continuous 
function f{r,t) of space r and time t satisfies the partial 
differential equation: 



+ wl{-iVr)f{z,t)=0, 



(6) 



and a;^(— «Vr) is the operator obtained after the sub- 
stitution k — > — iVr in the Fourier expansion of w^(k). 

It should be noted that the differential equation for 
the continuous envelope f{v,t) of the Wannier function 



wavepacket [Eq.(4)], as given by Eq.(6), is exact, and for 
any band of the PC an envelope equation can be writ- 
ten, the specific details of the band entering both in the 
dispersion curve w^(k) and in the shape of the corre- 
sponding Wannier function W„ [Eq.(2)]. 



III. SPATIAL AND SPATIO-TEMPORAL 
LOCALIZED WAVES 

The most general solution to the Wannier-function en- 
velope equation (6) is given by a superposition of func- 
tions ip{r,±t), where tp{r,t) is a solution to the wave 
equation: 



.dip , _ . , 

I— = UJn{-lVr)W- 



(7) 



We are know interested on the search for localized so- 
lutions to Eq.(7) such that corresponds to a wave 
propagating undistortcd with a group velocity Vg. To 
this aim, let us set -(/"(r, t) — g(r, t) exp(iko — i^lt), where 
ko is chosen inside the first Brillouin zone in the recipro- 
cal space and the frequency fl is chosen close to (but not 
necessarily coincident with) loq = a;„(ko). The envelope 
g then satisfies the wave equation 



.dg 



I— = [a;„(ko - iVr) - g. 



(8) 



We first note that, if g varies slowly with respect to 
the spatial variables r, at leading order one can ex- 
pand a;„(ko — iVr) up to first order around ko; tak- 
ing f2 = Wo, one obtains dg/dt + Vk^^'n • Vrg = 0, i.e. 
one retrieves the well-known result for which an arbitrary 
3D spatially-localized wavepacket travels undistorted, at 
leading order, with a group velocity given by Vk<^'n. Nev- 
ertheless, higher-order terms are generally responsible for 
wavepacket spreading, both in space and time. In order 
to find propagation-invariant envelope waves even when 
dispersive terms are accounted for, let us assume, without 
loss of generality, that {duin/dky)^^ = (9wn/9fcz)ko = 0, 
i.e. let us choose the orientation of the x axis such that 
the wavepacket group velocity VkWn is directed along 
this axis, and let us look for a propagation-invariant 
solution to Eq.(8) of the form g = g{x\,X2-iXz), with 
x\ = X — Vgt, X2 = y and xs = z, traveling along the x 
axis with a group velocity Vg , which is left undetermined 
at this stage. The function g then satisfies the following 
equation: 



-Wg— = [a;„(ko - iVx) - ^] g, 
whose solution can be written formally as: 



(9) 



gixi,X2,X3)^ y ci(32C?Q3 G(Q2,Q3)exp(iQ -x). (10) 

In Eq.(lO), X = (xi = x — Vgt,X2 = y,X3 = z), Q = 
(Qi, (52; Qa), G is an arbitrary spectral amplitude, and 
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Qi = Qi{Q2iQ?i) is implicitly defined by the following 
dispersion relation: 



w„(ko + Q) - fi - VgQi = 0. 



(11) 



To avoid the occurrence of evanescent (exponentially- 
growing) waves, the integral in Eq.(lO) is extended over 
the values of {Q2^ Qa) such that Qi, obtained after solv- 
ing Eq.(ll), turns out to be real-valued. We note that, 
for an arbitrary spectral amplitude G, Eq.(lO) represents 
an exact solution of the Wannier-function envelope equa- 
tion, which propagates undistorted with a group veloc- 
ity Vg, once the proper band dispersion curve w„(k) of 
the PC and corresponding dispersion relation (11) are 
computed, e.g. by numerical methods. For some spe- 
cific choices of the spectral amplitude G, in addition to 
undistorted wave propagation a certain degree of spatio- 
temporal wave localization can be obtained. It is worth 
to get some explicit examples, though approximate, of 
such 3D localized waves, admitting the integral repre- 
sentation given by Eq.(lO), and relate them to alreeidy 
known localized solutions to canonical wave equations }T\ . 
To this aim, we develop an asymptotic analysis of Eq.(ll) 
by assuming that the spectral amplitude G is nonvanish- 
ing in a narrow interval around Q2 = Q3 = 0, so that, for 
n close to Wo, the value of Qi, as obtained form Eq.(ll), 
is also close to Qi = 0. In this case, an approximate ex- 
pression for the dispersion relation Qi = Qi{Q2, Q3) can 
be obtained by expanding in Eq.(ll) the band dispersion 
curve w„(ko -I- Q) at around ko. We should distinguish 
two cases, depending on the value of the group velocity 
Vg, which is basically a free parameter in our analysis. 
First case. The first case corresponds to the choice of a 
group velocity Vg different from (and enough far form) 
dwn/dkx- In this case, the leading-order terms entering 
in Eq.(ll) after a power expansion of a;„(ko -I- Q) are 
quadratic in Q2, Q3 and linear in Qi; precisely, one has: 



If - Q.+^o-n+l ± i^Q.Q, - 0, (12) 



where fc, 



2 dkidkj 
for i = 1,2,3 and the derivatives of 



the band dispersion curve are calculated at k = kg. If 
the approximate expression of Qi, given Eq.(12), is intro- 
duced into Eq.(lO), one can easily show that the envelope 
g{xi,X2,X3) satisfies the difi^erential equation: 



afcT 



dxi 2 ^ dhdk, dx,dx 

i,j=2 





(13) 

Since the matrix d^tOn/dkidkj is symmetric, after a suit- 
able rotation of the {x2,X3) axes by the transformation 
Xj = TZjiXi (i,j = 2,3), where TZji is the orthogonal 
matrix that diagonalizes d^ujn/dkidkj, assuming with- 
out loss of generality Q — luq, Eq.(13) can be written in 
the canonical Schrodinger-like form: 

■ f _ ^ ^ 

\dk, dx. 



1 d'g 



1 d^g 

2"'ax;3' 



(14) 



where 012 and 013 are the eigenvalues of the 2x2 ma- 
trix d^LUn/dkidkj {i,j = 2,3). 3D localized waves to 
Eq.(14) are expressed in terms of well-known Gauss- 
Hermite functions, which are in general anisotropic for 
ck2 7^ ct3- These 3D localized waves, which exist re- 
gardless of the sign of 0:2 and a^, represent Gaussian- 
like beams, with exponential localization in the trans- 
verse (y, z) plane and algebraic localization, determined 
by the beam Rayleigh range, in the longitudinal x di- 
rection (and hence in time). These beams propagate 
undistorted along the x direction with an arbitrary group 
velocity Vg, either subluminal or superluminal, provided 
that Vg dujn/dkx- Such pulsed propagating Gaussian 
beams represent an extension, in a PC structure, of sim- 
ilar solutions found in vacuum (see |17| and references 
therein). In particular, the special case Vg = leads to 
stationary (monochromatic) Gaussian-like beams; note 
that the condition Vg ^ duon/dkx implies that such 
steady Gaussian beams do not exist in a PC close to a 
bandgap edge, where duJn/dkx vanishes. Other solutions 
to Eq.(14), leading to spatial 2D localized and monochro- 
matic waves in the transverse {y,z) plane (but delocalized 
in the longitudinal x direction) , can be search in the form 
g{xi,X2,X3) = s(a;2, X3) exp(iAa;i), where A is a propa- 
gation constant. If a2 and as have the same sign, the 
fimction s{x2,X3) satisfies a 2D Helmholtz equation, ad- 
mitting well-known Bessel-beam solutions in cylindrical 
coordinates. For a2 ^0^3, such solutions are anisotropic, 
and again they represent a generalization to a PC of well- 
known spatial Bessel beams in vacuum. If 012 and have 
opposite sign, one obtains a hyperbolic 2D equation (or, 
equivalently, a ID Klein-Gordon equation), which admits 
of 2D X-type localized solutions involving modified Bessel 
functions recently studied in Q (see Eqs.(3a) and (4) of 
Ref. 9]; see also 

Second case. The second case corresponds to the choice 
Vg = duJn/dkx. In this case, the leading-order approxi- 
mation to the dispersion relation [Eq.(ll)] should include 
also second-order derivatives with respect to xi of the 
band dispersion curve Ci;„(ko + Q), yielding: 



Wo 



1 ^ 



dkidkj 
1 



(15) 



where the derivatives of the band dispersion curve are 
calculated at k = ko. If the approximate expression of 
Qi, implicitly defined by the quadratic equation (15), is 
introduced into Eq.(lO), one can easily show that the 
envelope g(xi,X2,X3) satisfies this time the differential 
equation: 



{ujQ - n)g 



1 ^ 



^ dkidkj dxidxj 



(16) 



Since the matrix d^ujn/ dkidkj is symmetric, after a suit- 
able rotation of the {xi, X2, X3) axes by the transforma- 
tion Xj — TZjiXi {i,j = 1,2,3), where TZji is the orthogo- 
nal matrix that diagonalizes d^usn/dkidkj, Eq.(16) takes 
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the canonical form: 




where ai {i = 1, 2, 3) are the eigenvalues of the 3x3 ma- 
trix d'^uJn/dkidkj = 1,2,3). The sign of the eigen- 
values ai basically determines the elliptic or hyperbolic 
character of Eq.(17), and hence the nature of their solu- 
tions (see, e.g., [3). ai have the same sign, e.g. they 
are positive, for fl < ujq Eq.(17) reduces, after a scal- 
ing of axis length, to a 3D Hclmholtz equation, which in 
spherical coordinates admits of localized solutions in the 
form of sinc-shaped waves (see, e.g., 0,01). If, conversely, 
there is a sign discordance among the eigenvalues , one 
obtains a 2D Klein-Gordon equation, which admits of 
3D localized X-type waves which have been lengthly dis- 
cussed in many works (see, e.g., H, |E IS and references 
therein). In some special cases, one of the eigenvalues 
ai may vanish, which may yield further nonspreading 
wavepacket solutions. Notably, if ai = 0, the solution 
to Eq.(17) is given by 3(2:1,2:2,0:3) = h{xi)ip{x2, X3), 
where h is an arbitrary function of a;i = x — Vgt and 
ip satisfies a 2D Helmoltz equation for 02^3 > 0, ad- 
mitting Bessel beam solutions, or a ID Klein-Gordon 
equation for 020:3 < 0, admitting 2D X-type solutions. 
For these special solutions a cancellation of temporal dis- 
persion is attained. As the former case (020:3 > 0) ex- 
tends to a PC structure the so-called pulsed Bessel beams 
found in homogeneous dispersive media Q, the latter 
case (02O3 < 0) is rather peculiar for a PC structure, 
which realizes a bi-diffractive propagation regime 0, i.e. 



positive and negative diffraction along the two transverse 
directions y and z. Instead of pulses with a transverse 
Bessel beam profile, in this case one obtains a transverse 
X-shaped beam with an arbitrary longitudinal (tempo- 
ral) profile that propagates without spreading. 
As a final remark, we note that, though our analysis has 
been focused to a 3D PC, similar results can be obtained 
mutatis mutandis for the lower-dimensional case of a 2D 
PC. In this case, not considering out-of-plane propaga- 
tion, the fields depend solely on the two spatial variables 
X and y defining the PC plane, and Eqs.(14) and (17) are 
still valid provided that the terms involving the deriva- 
tives with respect to the 2:3 = z coordinate are dropped. 
In this case, Eq.(14) corresponds to a ID Schrodinger 
equation, whereas Eq.(17) corresponds to either a 2D 
Helmholtz equation or to a ID Klein-Gordon equation. 



IV. CONCLUSIONS 

In conclusion, a general analysis of wavepacket propa- 
gation in PCs, based on a Wannier function expansion 
approach, has been presented, and an exact envelope 
equation describing undistorted propagation of spatio- 
temporal localized waves has been derived. An asymp- 
totic analysis of the envelope equation shows that a wide 
class of localized (either spatial or spatio-temporal) waves 
exist, including propagating Gaussian beams, 2D and 3D 
X-type waves, sinc-shaped waves, pulsed Bessel beams 
and pulsed 2D X waves, some of which have been recently 
studied with reference to some specific models d, . 
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